function bdr=draw_matn(Y,X,sige,vecprM,prP)
% ==================================================
% function draw_matn
% Sample from the Matrivariate Normal distribution
% for the regression model
% Y=XB+E Var(E)=sige
% where a priori B~MN(vec(prM),prP^-1)
% and sige is (a priori) diffuse 
% Note that this implied independence 
% See for instance Kadiyala and Karlsson p. 106 and p. 111
% Input
% Y      [T nn] 
% X      [T nk] 
% sige   [nn nn] 
% vecprM [nn*nk 1]     Vec of the prior Mean 
% prP    [nn*nk nn*nk] Prior Precision 
%
% if vecprM and prP empty then do OLS regression
% Output 
% bdr    [1 nn*nk] draw 
% Alejandro Justiniano Sep 30 2008 
% =================================================
nk=size(X,2); nn=size(Y,2);
if nargin < 4 || isempty(vecprM)
    QQ=(X'*X)\eye(nk);
    bc=QQ*(X'*Y);
    % Regressors in columns for each Yf
    bdr=mvnrnd((bc(:))',kron(sige,QQ));
else
    % Posterior Variance
    % Recall prP is the prio precision
    invsigeXX=kron( (sige\eye(nn)) , (X'*X) );
    % Very large matrix inversion. In general cannot break into blocks 
    postV=(prP+invsigeXX)\eye(nk*nn);
    postV=0.5*(postV+postV'); 
    % OLS estimate
    bc=(X'*X)\(X'*Y);
    % Posterior Mean
    postM=postV*(prP*vecprM+invsigeXX*bc(:));
    bdr=mvnrnd((postM(:))',postV);
end